######################################################################################
# Project:  Complexity Project
#
# Task:     Generate fixed and random effects models
#
# Author:   Martijn Schoonvelde
# Date:     November 2021
#
######################################################################################

###############
#load libraries
###############

rm(list=ls())
library(lme4)
library(ggplot2)

#setwd(str_split_fixed(getwd(),"/Rep",2)[,1])
load("speeches.RData")

## Fixed-effect models
vars3 <- as.formula(complexity ~ lr + pc + year2 + office + factor(speaker))
reg.output.fixed <- list()
reg.output.fixed[[1]]  <- lm(vars3, data=speeches, subset=country=="DK" & data.set=="congress speeches")
reg.output.fixed[[2]]  <- lm(vars3, data=speeches, subset=country=="NL" & data.set=="congress speeches")

reg.output.fixed[[3]] <- lm(vars3, data=speeches, subset=country=="Germany" & data.set=="parliamentary speeches")
reg.output.fixed[[4]] <- lm(vars3, data=speeches, subset=country=="Spain"  & data.set=="parliamentary speeches") 
reg.output.fixed[[5]] <- lm(vars3, data=speeches, subset=country=="UK" & data.set=="parliamentary speeches")
reg.output.fixed[[6]] <- lm(vars3, data=speeches, subset=country=="Sweden" & data.set=="parliamentary speeches")
reg.output.fixed[[7]] <-lm(vars3, data=speeches, subset=country=="Netherlands" & data.set=="parliamentary speeches") 

#plot data
##########
coefs <- unlist(lapply(reg.output.fixed, function(x) coef(x)[3]))
se <- unlist(lapply(reg.output.fixed, function(x) sqrt(diag(vcov(x)))[3]))
min <- coefs - 1.96*se
max <- coefs + 1.96*se
names <- c("Party congresses DNK", "Party congresses NLD", "Parliament GER", "Parliament ESP", "Parliament GBR", "Parliament SWE","Parliament NLD") 
plot.data <- data.frame(names,coefs,se,min,max)
order.names <- c(2,1,7,6,5,4,3)
plot.data <- plot.data[order.names,]
plot.data$names <- factor(plot.data$names, levels=plot.data$names)
plot.data$institution <- c(rep("Congress","2"), rep("Parliament","5"))
plot.data$institution <- factor(plot.data$institution, levels = c("Parliament", "Congress"))
plot.data$model <- "fixed effects"

coef.plot3 <- ggplot(plot.data, aes(x=names, y=coefs)) +
  geom_point(size = 0.75) +
  geom_errorbar(ymin=plot.data$min[],ymax=plot.data$max[], width=.15) +
  coord_flip() + 
  ylab("Fixed effect of conservatism") +
  xlab("") + 
  geom_hline(yintercept=0) + theme(panel.spacing = unit(2, "lines")) +  ylim(c(-0.45,0.40)) +
  theme_minimal() + facet_wrap(institution~., scales="free",nrow = 2) 


## Random-intercept models
vars3 <- as.formula(complexity ~ lr + pc + year2 + office + (1 | speaker))
reg.output.random <- list()
reg.output.random[[1]]  <- lmer(vars3, data=speeches, subset=country=="DK" & data.set=="congress speeches")
reg.output.random[[2]]  <- lmer(vars3, data=speeches, subset=country=="NL" & data.set=="congress speeches")

reg.output.random[[3]] <- lmer(vars3, data=speeches, subset=country=="Germany" & data.set=="parliamentary speeches")
reg.output.random[[4]] <- lmer(vars3, data=speeches, subset=country=="Spain"  & data.set=="parliamentary speeches") 
reg.output.random[[5]] <- lmer(vars3, data=speeches, subset=country=="UK" & data.set=="parliamentary speeches")
reg.output.random[[6]] <- lmer(vars3, data=speeches, subset=country=="Sweden" & data.set=="parliamentary speeches")
reg.output.random[[7]] <-lmer(vars3, data=speeches, subset=country=="Netherlands" & data.set=="parliamentary speeches") 


coef_1 <- coef(reg.output.random[[1]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[1]]))
min_1 <- confint[5,1]
max_1 <- confint[5,2]

coef_2 <- coef(reg.output.random[[2]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[2]]))
min_2 <- confint[5,1]
max_2 <- confint[5,2]


coef_3 <- coef(reg.output.random[[3]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[3]]))
min_3 <- confint[5,1]
max_3 <- confint[5,2]


coef_4 <- coef(reg.output.random[[4]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[4]]))
min_4 <- confint[5,1]
max_4 <- confint[5,2]

coef_5 <- coef(reg.output.random[[5]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[5]]))
min_5 <- confint[5,1]
max_5 <- confint[5,2]

coef_6 <- coef(reg.output.random[[6]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[6]]))
min_6 <- confint[5,1]
max_6 <- confint[5,2]

coef_7 <- coef(reg.output.random[[7]])$speaker$pc[1]
confint <- as.data.frame(confint(reg.output.random[[7]]))
min_7 <- confint[5,1]
max_7 <- confint[5,2]


coefs <- c(coef_1, coef_2, coef_3, coef_4, coef_5, coef_6, coef_7)
min <- c(min_1, min_2, min_3, min_4, min_5, min_6, min_7)
max <- c(max_1, max_2, max_3, max_4, max_5, max_6, max_7)


names <- c("Party congresses DNK", "Party congresses NLD", "Parliament GER", "Parliament ESP", "Parliament GBR", "Parliament SWE","Parliament NLD") 
plot.data.re <- data.frame(names,coefs,se,min,max)
plot.data.re$se <- NA
order.names <- c(2,1,7,6,5,4,3)
plot.data.re <- plot.data.re[order.names,]
plot.data.re$names <- factor(plot.data.re$names, levels=plot.data.re$names)
plot.data.re$institution <- c(rep("Congress","2"), rep("Parliament","5"))
plot.data.re$institution <- factor(plot.data.re$institution, levels = c("Parliament", "Congress"))

plot.data.re$model <- "random effects"

plot.data.total <- rbind(plot.data, plot.data.re)

coef.plot3 <- ggplot(plot.data.re, aes(x=names, y=coefs)) +
  geom_point(size = 0.75) +
  geom_errorbar(ymin=plot.data$min[],ymax=plot.data$max[], width=.15) +
  coord_flip() + 
  ylab("Fixed effect of conservatism") +
  xlab("") + 
  geom_hline(yintercept=0) + theme(panel.spacing = unit(2, "lines")) +  ylim(c(-0.45,0.40)) +
    theme_minimal() + facet_wrap(institution~., scales="free",nrow = 2)


total.sample.plot <- ggplot(plot.data.total[], aes(x=names, 
                                            y=coefs, 
                                            ymin = min,
                                            ymax = max)) +
  geom_pointrange(position = position_dodge(width = 0.5)) +
  coord_flip() + facet_wrap(~ model) +
  ylab("Effect of lib-cons ideology on F-K") +
  xlab("Corpora") + 
 # scale_y_continuous(limits = c(-0.45, 0.25)) + 
  geom_hline(yintercept=0) +
  theme_minimal()


ggsave(total.sample.plot , file = "fe_and_re.png",
       width = 10, height = 8)
